#ifndef VtkWriter_H
#define VtkWriter_H

#include "pugixml/pugixml.hpp"
#include "pugixml/pugiconfig.hpp"

#include "Mesh/Mesh.h"
#include "Mesh/ParallelMesh.h"
#include "Config/Config.h"
#include "IO/IO.h"

#include <string>
#include <vector>
#include <csignal>

using namespace std;
using namespace pugi;

class VtkWriter
{
public:

    VtkWriter(Config* configInp, IO* ioInp)
    {
        config = configInp;
        mesh = config->mesh;
        io = ioInp;

        globalMesh = mesh->serialMesh;
    };

    ~VtkWriter();

    void writeVtuMesh()
    {
        std::string vtuFileName = "result/vtk/mesh.vtu." + std::to_string(config->pmpi->get_rank());
        if (config->pmpi->get_size() == 1)
        {
            vtuFileName = "result/vtk/mesh.vtu";
        }

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto unstructuredGridNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = unstructuredGridNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(mesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(mesh->volumeList.size()).c_str();

        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = mesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());


        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            EasyMesh::Volume& volume = mesh->volumeList[iCell];
            connectivityDataString += to_string(volume.pointList[0]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[1]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[2]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[3]);
            connectivityDataString += "\n";
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            iOffset += 4;
            offsetsDataString += to_string(iOffset);
            offsetsDataString += "\n";
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "UInt8";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";
        string typesDataString = "\n";
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            typesDataString += "10";
            typesDataString += "\n";
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());

        auto pointDataNode = pieceNode.append_child("PointData");

        auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
        pointDataNodeDataArrayPotential.append_attribute("Name") = "Potential";
        pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
        pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
        pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
        string potentialDataString = "\n";
        for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = mesh->pointList[iPoint];
            //potentialDataString += to_string(iPoint);
            potentialDataString += to_string(point.position[0]);
            potentialDataString += "\n";
        }
        pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());


        auto cellDataNode = pieceNode.append_child("CellData");

        doc.save_file(vtuFileName.c_str());

    };

    void writeVtuField(string timeString)
    {
        if (config->dimension == 3)
        {
            writeVtuField3d(timeString);
        }
        else if (config->dimension == 2)
        {
            writeVtuField2d(timeString);
        }
        else
        {
            log_error << "the dimension not support: " << config->dimension;
            exit(0);
        }
    }

    void writeVtuField3d(string timeString)
    {
        std::string vtuFileName = "result/vtk/fieldProc" + std::to_string(config->pmpi->get_rank()) + "_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto unstructuredGridNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = unstructuredGridNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(mesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(mesh->volumeList.size()).c_str();

        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = mesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());


        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            EasyMesh::Volume& volume = mesh->volumeList[iCell];
            connectivityDataString += to_string(volume.pointList[0]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[1]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[2]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[3]);
            connectivityDataString += "\n";
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            iOffset += 4;
            offsetsDataString += to_string(iOffset);
            offsetsDataString += "\n";
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "UInt8";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";
        string typesDataString = "\n";
        for (int iCell = 0; iCell < mesh->volumeList.size(); iCell++)
        {
            typesDataString += "10";
            typesDataString += "\n";
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());


        auto pointDataNode = pieceNode.append_child("PointData");

        //write scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = (io->fields[iField]->name + "_oneStep").c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                //potentialDataString += to_string(iPoint);
                potentialDataString += to_string(io->fields[iField]->data[iPoint]);
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }

        //write average scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                //potentialDataString += to_string(iPoint);
                potentialDataString += to_string(io->fields[iField]->data_average[iPoint]);
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }

        //write vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = (io->vector_fields[iField]->name + "_oneStep").c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }

        //write average vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->vector_fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }



        //write cell data
        auto cellDataNode = pieceNode.append_child("CellData");

        doc.save_file(vtuFileName.c_str());

    };

    void writeVtuField2d(string timeString)
    {
        std::string vtuFileName = "result/vtk/fieldProc" + std::to_string(config->pmpi->get_rank()) + "_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto unstructuredGridNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = unstructuredGridNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(mesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(mesh->numberCell).c_str();

        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = mesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());


        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (int iCell = 0; iCell < mesh->numberCell; iCell++)
        {
            EasyMesh::Face& face = mesh->faceList[iCell];
            connectivityDataString += to_string(face.pointList[0]);
            connectivityDataString += " ";
            connectivityDataString += to_string(face.pointList[1]);
            connectivityDataString += " ";
            connectivityDataString += to_string(face.pointList[2]);
            connectivityDataString += "\n";
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (int iCell = 0; iCell < mesh->numberCell; iCell++)
        {
            iOffset += 3;
            offsetsDataString += to_string(iOffset);
            offsetsDataString += "\n";
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "UInt8";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";
        string typesDataString = "\n";
        for (int iCell = 0; iCell < mesh->numberCell; iCell++)
        {
            typesDataString += "5";
            typesDataString += "\n";
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());


        auto pointDataNode = pieceNode.append_child("PointData");

        //write scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                //potentialDataString += to_string(iPoint);
                potentialDataString += to_string(io->fields[iField]->data[iPoint]);
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }

        //write average scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = (io->fields[iField]->name + "Average").c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                //potentialDataString += to_string(iPoint);
                potentialDataString += to_string(io->fields[iField]->data_average[iPoint]);
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }

        //write vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->vector_fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }

        //write average vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = (io->vector_fields[iField]->name + "Average").c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }



        //write cell data
        auto cellDataNode = pieceNode.append_child("CellData");

        doc.save_file(vtuFileName.c_str());

    };



    //write average global field for visulization

    void writeVtuFieldGlobal(string timeString)
    {
        if (config->dimension == 3)
        {
            writeVtuFieldGlobal3d(timeString);
        }
        else if (config->dimension == 2)
        {
            writeVtuFieldGlobal2d(timeString);
        }
        else
        {
            log_error << "the dimension not support: " << config->dimension;
            exit(0);
        }
    }


    void writeVtuFieldGlobal3d(string timeString)
    {
        std::string vtuFileName = "result/vtk/field_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto unstructuredGridNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = unstructuredGridNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(globalMesh->volumeList.size()).c_str();

        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());


        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (int iCell = 0; iCell < globalMesh->volumeList.size(); iCell++)
        {
            EasyMesh::Volume& volume = globalMesh->volumeList[iCell];
            connectivityDataString += to_string(volume.pointList[0]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[1]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[2]);
            connectivityDataString += " ";
            connectivityDataString += to_string(volume.pointList[3]);
            connectivityDataString += "\n";
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (int iCell = 0; iCell < globalMesh->volumeList.size(); iCell++)
        {
            iOffset += 4;
            offsetsDataString += to_string(iOffset);
            offsetsDataString += "\n";
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "UInt8";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";
        string typesDataString = "\n";
        for (int iCell = 0; iCell < globalMesh->volumeList.size(); iCell++)
        {
            typesDataString += "10";
            typesDataString += "\n";
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());


        auto pointDataNode = pieceNode.append_child("PointData");


        //write global scalar data
        if (config->debugLevel > 1)
        {
            for (int iField = 0; iField < io->fields.size(); iField++)
            {
                auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
                pointDataNodeDataArrayPotential.append_attribute("Name") = (io->fields[iField]->name + "_OneStep").c_str();
                pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
                pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
                pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
                string potentialDataString = "\n";
                for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
                {
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << io->fields[iField]->data_global[iPoint];
                    ss >> str;
                    potentialDataString += str;
                    potentialDataString += "\n";
                }
                pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
            }
        }



        //write averge global scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = (io->fields[iField]->name).c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
            {
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->fields[iField]->data_average_global[iPoint];
                ss >> str;
                potentialDataString += str;
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }


        //write global vector data
        if (config->debugLevel > 1)
        {
            for (int iField = 0; iField < io->vector_fields.size(); iField++)
            {
                auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
                pointDataNodeDataArrayPotential.append_attribute("Name") = (io->vector_fields[iField]->name + "_OneStep").c_str();
                pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
                pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
                pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
                string pointDataString = "\n";
                for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
                {
                    EasyMesh::Point& point = globalMesh->pointList[iPoint];
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[0];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[1];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[2];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += "\n";
                }
                pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
            }
        }
        


        //write average global vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->vector_fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = globalMesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }





        //write cell data
        auto cellDataNode = pieceNode.append_child("CellData");

        doc.save_file(vtuFileName.c_str());

    };

    void writeVtuFieldGlobal2d(string timeString)
    {
        std::string vtuFileName = "result/vtk/field_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto unstructuredGridNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = unstructuredGridNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(globalMesh->faceList.size()).c_str();

        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());


        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (int iCell = 0; iCell < globalMesh->numberCell; iCell++)
        {
            EasyMesh::Face& face = globalMesh->faceList[iCell];
            connectivityDataString += to_string(face.pointList[0]);
            connectivityDataString += " ";
            connectivityDataString += to_string(face.pointList[1]);
            connectivityDataString += " ";
            connectivityDataString += to_string(face.pointList[2]);
            connectivityDataString += "\n";
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (int iCell = 0; iCell < globalMesh->numberCell; iCell++)
        {
            iOffset += 3;
            offsetsDataString += to_string(iOffset);
            offsetsDataString += "\n";
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "UInt8";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";
        string typesDataString = "\n";
        for (int iCell = 0; iCell < globalMesh->numberCell; iCell++)
        {
            typesDataString += "5";
            typesDataString += "\n";
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());


        auto pointDataNode = pieceNode.append_child("PointData");


        //write global scalar data
        if (config->debugLevel > 1)
        {
            for (int iField = 0; iField < io->fields.size(); iField++)
            {
                auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
                pointDataNodeDataArrayPotential.append_attribute("Name") = (io->fields[iField]->name + "_oneStep").c_str();
                pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
                pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
                pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
                string potentialDataString = "\n";
                for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
                {
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << io->fields[iField]->data_global[iPoint];
                    ss >> str;
                    potentialDataString += str;
                    potentialDataString += "\n";
                }
                pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
            }
        }


        //write averge global scalar data
        for (int iField = 0; iField < io->fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "1";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string potentialDataString = "\n";
            for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
            {
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->fields[iField]->data_average_global[iPoint];
                ss >> str;
                potentialDataString += str;
                potentialDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(potentialDataString.c_str());
        }

        //write global vector data
        if (config->debugLevel > 1)
        {
            for (int iField = 0; iField < io->vector_fields.size(); iField++)
            {
                auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
                pointDataNodeDataArrayPotential.append_attribute("Name") = (io->vector_fields[iField]->name + "_oneStep").c_str();
                pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
                pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
                pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
                string pointDataString = "\n";
                for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
                {
                    EasyMesh::Point& point = globalMesh->pointList[iPoint];
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[0];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[1];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << io->vector_fields[iField]->data_global[iPoint].x[2];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += "\n";
                }
                pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
            }
        }
        



        //write average global vector data
        for (int iField = 0; iField < io->vector_fields.size(); iField++)
        {
            auto pointDataNodeDataArrayPotential = pointDataNode.append_child("DataArray");
            pointDataNodeDataArrayPotential.append_attribute("Name") = io->vector_fields[iField]->name.c_str();
            pointDataNodeDataArrayPotential.append_attribute("type") = "Float64";
            pointDataNodeDataArrayPotential.append_attribute("NumberOfComponents") = "3";
            pointDataNodeDataArrayPotential.append_attribute("format") = "ascii";
            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = globalMesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << io->vector_fields[iField]->data_average_global[iPoint].x[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointDataNodeDataArrayPotential.text().set(pointDataString.c_str());
        }





        //write cell data
        auto cellDataNode = pieceNode.append_child("CellData");

        doc.save_file(vtuFileName.c_str());

    };


    void writeVtpMesh()
    {
        if (config->dimension == 2)
        {
            return;
        }

        //write faceGroupList of mesh, one faceGroup(boundary face) to one vtp file

        for (auto iter = mesh->faceGroupMap.begin(); iter != mesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;


            std::string vtpFileName = "result/vtk/meshFaceGroup_" + faceGroupName + ".vtp." + std::to_string(config->pmpi->get_rank());
            if (config->pmpi->get_size() == 1)
            {
                vtpFileName = "result/vtk/meshFaceGroup_" + faceGroupName + ".vtp";
            }

            xml_document doc;
            auto vtkFileNode = doc.append_child("VTKFile");
            vtkFileNode.append_attribute("type") = "PolyData";

            auto polyDataNode = vtkFileNode.append_child("PolyData");

            auto pieceNode = polyDataNode.append_child("Piece");
            pieceNode.append_attribute("NumberOfPoints") = to_string(mesh->pointList.size()).c_str();
            pieceNode.append_attribute("NumberOfVerts") = "0";
            pieceNode.append_attribute("NumberOfLines") = "0";
            pieceNode.append_attribute("NumberOfStrips") = "0";
            pieceNode.append_attribute("NumberOfPolys") = to_string(faceGroup.faceList.size()).c_str();

            auto pointsNode = pieceNode.append_child("Points");
            auto pointsNodeDataArray = pointsNode.append_child("DataArray");
            pointsNodeDataArray.append_attribute("type") = "Float64";
            pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
            pointsNodeDataArray.append_attribute("format") = "ascii";

            string pointDataString = "\n";
            for (int iPoint = 0; iPoint < mesh->pointList.size(); iPoint++)
            {
                EasyMesh::Point& point = mesh->pointList[iPoint];
                std::stringstream ss;
                std::string str;

                ss.clear();
                ss << std::setprecision(14) << point.position[0];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << point.position[1];
                ss >> str;
                pointDataString += str;
                pointDataString += " ";

                ss.clear();
                ss << std::setprecision(14) << point.position[2];
                ss >> str;
                pointDataString += str;
                pointDataString += "\n";
            }
            pointsNodeDataArray.text().set(pointDataString.c_str());


            auto polysNode = pieceNode.append_child("Polys");

            auto polysNodeDataArrayConnectivity = polysNode.append_child("DataArray");
            polysNodeDataArrayConnectivity.append_attribute("type") = "Int32";
            polysNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
            polysNodeDataArrayConnectivity.append_attribute("format") = "ascii";

            string connectivityDataString = "\n";
            for (int i = 0; i < faceGroup.faceList.size(); i++)
            {
                int iFace = faceGroup.faceList[i];
                EasyMesh::Face& face = mesh->faceList[iFace];
                connectivityDataString += to_string(face.pointList[0]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[1]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[2]);
                connectivityDataString += "\n";
            }
            polysNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

            auto polysNodeDataArrayOffsets = polysNode.append_child("DataArray");
            polysNodeDataArrayOffsets.append_attribute("type") = "Int32";
            polysNodeDataArrayOffsets.append_attribute("Name") = "offsets";
            polysNodeDataArrayOffsets.append_attribute("format") = "ascii";

            string offsetsDataString = "\n";
            int iOffset = 0;
            for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
            {
                iOffset += 3;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
            polysNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

            //PolyData存储的是面，以及面上的数据，因此这里的cell指的是Polys表示的面，值在面心
            auto cellDataNode = pieceNode.append_child("CellData");

            auto cellDataNodeNodeDataArrayParticleFlux = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayParticleFlux.append_attribute("Name") = "ParticleFlux";
            cellDataNodeNodeDataArrayParticleFlux.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayParticleFlux.append_attribute("NumberOfComponents") = "1";
            cellDataNodeNodeDataArrayParticleFlux.append_attribute("format") = "ascii";
            string particleFluxDataString = "\n";
            for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
            {
                //Point& point = mesh->pointList[iPoint];
                //potentialDataString += to_string(iPoint);
                particleFluxDataString += to_string(iFace);
                particleFluxDataString += "\n";
            }
            cellDataNodeNodeDataArrayParticleFlux.text().set(particleFluxDataString.c_str());

            doc.save_file(vtpFileName.c_str());

        }



    };


    void writeVtpSpecies(string timeString)
    {
        //write species particle position velocity to one vtp file

        for (int i_species = 0; i_species < io->speciesCollection->species_list.size(); i_species++)
        {
            auto species = io->speciesCollection->species_list[i_species];
            int particleNumber = species->get_particle_number();

            //log_info << "begin writing species333 " << species->species_param.isWriteVtk;
            if (!species->species_param.isWrite)
            {
                continue;
            }

            std::string vtpFileName = "result/vtk/species_" + species->species_param.species_name + "Proc" + std::to_string(config->pmpi->get_rank()) + "_" + timeString + ".vtp";
            if (config->pmpi->get_size() == 1)
            {
                vtpFileName = "result/vtk/species_" + species->species_param.species_name + "_" + timeString + ".vtp";
            }


            log_info << "writing species vtk file: " << species->species_param.species_name;

            xml_document doc;
            auto vtkFileNode = doc.append_child("VTKFile");
            vtkFileNode.append_attribute("type") = "PolyData";

            auto polyDataNode = vtkFileNode.append_child("PolyData");

            auto pieceNode = polyDataNode.append_child("Piece");
            pieceNode.append_attribute("NumberOfPoints") = to_string(species->get_particle_number()).c_str();
            pieceNode.append_attribute("NumberOfVerts") = to_string(species->get_particle_number()).c_str();
            pieceNode.append_attribute("NumberOfLines") = "0";
            pieceNode.append_attribute("NumberOfStrips") = "0";
            pieceNode.append_attribute("NumberOfPolys") = "0";

            auto pointsNode = pieceNode.append_child("Points");
            auto pointsNodeDataArray = pointsNode.append_child("DataArray");
            pointsNodeDataArray.append_attribute("type") = "Float64";
            pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
            pointsNodeDataArray.append_attribute("format") = "ascii";

            //log_info << "writing species position ";
            string pointDataString = "\n";
            for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
            {
                auto& cell_particles = species->particles[i_volume];
                for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                {
                    auto& particle = cell_particles.data[i_particle];
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[0];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[1];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[2];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += "\n";
                }
            }
            pointsNodeDataArray.text().set(pointDataString.c_str());
            //log_info << "writing species position 1111";

            auto vertsNode = pieceNode.append_child("Verts");

            auto vertsNodeDataArrayConnectivity = vertsNode.append_child("DataArray");
            vertsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
            vertsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
            vertsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

            //log_info << "writing connectivity 1111";
            string connectivityDataString = "\n";
            for (int i = 0; i < particleNumber; i++)
            {
                //log_info << "writing species position 2222  " <<i;
                connectivityDataString += to_string(i);
                connectivityDataString += "\n";
            }
            vertsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());
            //log_info << "writing connectivity 2222";

            auto vertsNodeDataArrayOffsets = vertsNode.append_child("DataArray");
            vertsNodeDataArrayOffsets.append_attribute("type") = "Int32";
            vertsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
            vertsNodeDataArrayOffsets.append_attribute("format") = "ascii";

            string offsetsDataString = "\n";
            int iOffset = 0;
            for (int i = 0; i < particleNumber; i++)
            {
                iOffset += 1;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
            vertsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

            //velocity
            auto cellDataNode = pieceNode.append_child("CellData");

            //log_info << "writing species velocity ";
            auto cellDataNodeNodeDataArrayVelocity = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayVelocity.append_attribute("Name") = "velocity";
            cellDataNodeNodeDataArrayVelocity.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayVelocity.append_attribute("NumberOfComponents") = "3";
            cellDataNodeNodeDataArrayVelocity.append_attribute("format") = "ascii";
            string velocityDataString = "\n";
            for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
            {
                auto& cell_particles = species->particles[i_volume];
                for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                {
                    auto& particle = cell_particles.data[i_particle];

                    velocityDataString += to_string(particle.velocity[0]);
                    velocityDataString += " ";
                    velocityDataString += to_string(particle.velocity[1]);
                    velocityDataString += " ";
                    velocityDataString += to_string(particle.velocity[2]);
                    velocityDataString += "\n";
                }
            }
            cellDataNodeNodeDataArrayVelocity.text().set(velocityDataString.c_str());

            doc.save_file(vtpFileName.c_str());

            log_info << "done writing species: " << species->species_param.species_name;

        }



    };

    void writeVtuSpecies(string timeString)
    {
        //write species particle position velocity to multi vtu file
        for (int i_species = 0; i_species < io->speciesCollection->species_list.size(); i_species++)
        {
            auto species = io->speciesCollection->species_list[i_species];
            int particleNumber = species->get_particle_number();

            //log_info << "begin writing species333 " << species->species_param.isWriteVtk;
            if (!species->species_param.isWrite)
            {
                continue;
            }

            std::string vtpFileName = "result/vtk/species_" + species->species_param.species_name + "Proc" + std::to_string(config->pmpi->get_rank()) + "_" + timeString + ".vtu";
            if (config->pmpi->get_size() == 1)
            {
                vtpFileName = "result/vtk/species_" + species->species_param.species_name + "_" + timeString + ".vtu";
            }


            log_info << "writing species vtk file: " << species->species_param.species_name;

            xml_document doc;
            auto vtkFileNode = doc.append_child("VTKFile");
            vtkFileNode.append_attribute("type") = "UnstructuredGrid";

            auto polyDataNode = vtkFileNode.append_child("UnstructuredGrid");

            auto pieceNode = polyDataNode.append_child("Piece");
            pieceNode.append_attribute("NumberOfPoints") = to_string(species->get_particle_number()).c_str();
            pieceNode.append_attribute("NumberOfCells") = to_string(species->get_particle_number()).c_str();

            auto pointsNode = pieceNode.append_child("Points");
            auto pointsNodeDataArray = pointsNode.append_child("DataArray");
            pointsNodeDataArray.append_attribute("type") = "Float64";
            pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
            pointsNodeDataArray.append_attribute("format") = "ascii";

            //log_info << "writing species position ";
            string pointDataString = "\n";
            for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
            {
                //log_info << "species11 begin, " << i_volume;
                auto& cell_particles = species->particles[i_volume];
                for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                {
                    //log_info << "species begin, " << i_particle << " " << cell_particles.data.size();
                    auto& particle = cell_particles.data[i_particle];
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[0];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[1];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << particle.position[2];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += "\n";
                    //log_info << "species end, " << i_particle;
                }
                //log_info << "species11 end, " << i_volume;
            }
            pointsNodeDataArray.text().set(pointDataString.c_str());
            //log_info << "writing species position 1111";

            auto vertsNode = pieceNode.append_child("Cells");

            auto vertsNodeDataArrayConnectivity = vertsNode.append_child("DataArray");
            vertsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
            vertsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
            vertsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

            //log_info << "writing connectivity 1111";
            string connectivityDataString = "\n";
            for (int i = 0; i < particleNumber; i++)
            {
                //log_info << "writing species position 2222  " <<i;
                connectivityDataString += to_string(i);
                connectivityDataString += "\n";
            }
            vertsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());
            //log_info << "writing connectivity 2222";

            //offsets
            auto vertsNodeDataArrayOffsets = vertsNode.append_child("DataArray");
            vertsNodeDataArrayOffsets.append_attribute("type") = "Int32";
            vertsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
            vertsNodeDataArrayOffsets.append_attribute("format") = "ascii";

            string offsetsDataString = "\n";
            int iOffset = 0;
            for (int i = 0; i < particleNumber; i++)
            {
                iOffset += 1;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
            vertsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

            //Types
            auto vertsNodeDataArrayTypes = vertsNode.append_child("DataArray");
            vertsNodeDataArrayTypes.append_attribute("type") = "Int32";
            vertsNodeDataArrayTypes.append_attribute("Name") = "types";
            vertsNodeDataArrayTypes.append_attribute("format") = "ascii";

            string typesDataString = "\n";
            for (int i = 0; i < particleNumber; i++)
            {
                typesDataString += "1";
                typesDataString += "\n";
            }
            vertsNodeDataArrayTypes.text().set(typesDataString.c_str());



            //velocity
            auto cellDataNode = pieceNode.append_child("PointData");

            //log_info << "writing species velocity " << vtpFileName;
            auto cellDataNodeNodeDataArrayVelocity = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayVelocity.append_attribute("Name") = "velocity";
            cellDataNodeNodeDataArrayVelocity.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayVelocity.append_attribute("NumberOfComponents") = "3";
            cellDataNodeNodeDataArrayVelocity.append_attribute("format") = "ascii";
            string velocityDataString = "\n";
            for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
            {
                auto& cell_particles = species->particles[i_volume];
                for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                {
                    auto& particle = cell_particles.data[i_particle];

                    velocityDataString += to_string(particle.velocity[0]);
                    velocityDataString += " ";
                    velocityDataString += to_string(particle.velocity[1]);
                    velocityDataString += " ";
                    velocityDataString += to_string(particle.velocity[2]);
                    velocityDataString += "\n";
                }
            }
            cellDataNodeNodeDataArrayVelocity.text().set(velocityDataString.c_str());
            //log_info << "done writing species111: " << species->species_param.species_name << " " << vtpFileName;

            doc.save_file(vtpFileName.c_str());

            log_info << "done writing species: " << species->species_param.species_name;

        }



    };

    void writeVtuSpeciesToOneFile(string timeString)
    {
        //同步不同进程的粒子数，为了规约粒子到主进程
        std::vector<int> macroParticleNumber;
        macroParticleNumber.resize(io->speciesCollection->species_list.size() * config->pmpi->get_size(), 0);
        for (int iSpecies = 0; iSpecies < io->speciesCollection->species_list.size(); iSpecies++)
        {
            Species* species_ = io->speciesCollection->species_list[iSpecies];
            for (int iRank = 0; iRank < config->pmpi->get_size(); iRank++)
            {
                if (config->pmpi->get_rank() == iRank)
                {
                    macroParticleNumber[iSpecies * config->pmpi->get_size() + iRank] = species_->getMacroParticleNumber();
                }
                else
                {
                    macroParticleNumber[iSpecies * config->pmpi->get_size() + iRank] = 0.0;
                }
            }
        }
        std::vector<int> macroParticleNumberTemp = macroParticleNumber;
        config->pmpi->all_reduce_sum_int(macroParticleNumberTemp.data(), macroParticleNumber.data(), macroParticleNumber.size());

        std::vector<int> macroParticleNumberGlobal;
        macroParticleNumberGlobal.resize(io->speciesCollection->species_list.size());
        for (int i_species = 0; i_species < io->speciesCollection->species_list.size(); i_species++)
        {
            Species* species_ = io->speciesCollection->species_list[i_species];
            macroParticleNumberGlobal[i_species] = species_->getGlobalMacroParticleNumber();
        }



        //write species particle position velocity to one vtu file
        for (int i_species = 0; i_species < io->speciesCollection->species_list.size(); i_species++)
        {
            auto species = io->speciesCollection->species_list[i_species];

            int particleNumberGlobal = species->getGlobalMacroParticleNumber();
            int particleNumberLocal = species->getMacroParticleNumber();

            //log_info << "begin writing species333 " << species->species_param.isWriteVtk;
            if (!species->species_param.isWrite)
            {
                continue;
            }

            std::string vtpFileName = "result/vtk/species_" + species->species_param.species_name + "_" + timeString + ".vtu";


            if(config->debugLevel > 1) log_info << "writing species vtk file: " << species->species_param.species_name;


            xml_document doc;
            auto vtkFileNode = doc.append_child("VTKFile");
            vtkFileNode.append_attribute("type") = "UnstructuredGrid";

            auto polyDataNode = vtkFileNode.append_child("UnstructuredGrid");

            auto pieceNode = polyDataNode.append_child("Piece");
            pieceNode.append_attribute("NumberOfPoints") = to_string(particleNumberGlobal).c_str();
            pieceNode.append_attribute("NumberOfCells") = to_string(particleNumberGlobal).c_str();

            auto pointsNode = pieceNode.append_child("Points");
            auto pointsNodeDataArray = pointsNode.append_child("DataArray");
            pointsNodeDataArray.append_attribute("type") = "Float64";
            pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
            pointsNodeDataArray.append_attribute("format") = "ascii";

            //get all position from all processes
            if (config->debugLevel > 1) log_info << "writing species vtk file: gather postion";
            double* position = new double[particleNumberGlobal * 3];
            double* positionTemp = new double[particleNumberLocal * 3];
            for (int i = 0; i < particleNumberLocal * 3; i++)
            {
                positionTemp[i] = 0.0;
            }

            int* offset = new int[config->pmpi->get_size()];
            int* count = new int[config->pmpi->get_size()];
            offset[0] = 0;
            for (int iRank = 0; iRank < config->pmpi->get_size() - 1; iRank++)
            {
                offset[iRank+1] = offset[iRank] + 3 * macroParticleNumber[i_species * config->pmpi->get_size() + iRank];
            }
            for (int iRank = 0; iRank < config->pmpi->get_size(); iRank++)
            {
                count[iRank] = 3 * macroParticleNumber[i_species * config->pmpi->get_size() + iRank];
            }


            int iParticle = 0;
            for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
            {
                auto& cell_particles = species->particles[i_volume];
                for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                {
                    auto& particle = cell_particles.data[i_particle];
                    positionTemp[iParticle * 3 + 0] = particle.position[0];
                    positionTemp[iParticle * 3 + 1] = particle.position[1];
                    positionTemp[iParticle * 3 + 2] = particle.position[2];

                    iParticle++;
                }
            }

            config->pmpi->gatherv_double(positionTemp, 3*particleNumberLocal, position, count, offset, 0);


            //log_info << "writing species position ";
            if (config->debugLevel > 1) log_info << "writing species vtk file: write postion";
            if (config->pmpi->is_master())
            {
                string pointDataString = "\n";
                for (int iParticle = 0; iParticle < particleNumberGlobal; iParticle++)
                {
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << position[iParticle * 3 + 0];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << position[iParticle * 3 + 1];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += " ";

                    ss.clear();
                    ss << std::setprecision(14) << position[iParticle * 3 + 2];
                    ss >> str;
                    pointDataString += str;
                    pointDataString += "\n";
                    //log_info << "species end, " << i_particle;
                }
                pointsNodeDataArray.text().set(pointDataString.c_str());
                //log_info << "writing species position 1111";
            }
            delete[] position;
            delete[] positionTemp;
            delete[] count;
            delete[] offset;
            

            auto vertsNode = pieceNode.append_child("Cells");

            auto vertsNodeDataArrayConnectivity = vertsNode.append_child("DataArray");
            vertsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
            vertsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
            vertsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

            //log_info << "writing connectivity 1111";
            if (config->pmpi->is_master())
            {
                string connectivityDataString = "\n";
                for (int i = 0; i < particleNumberGlobal; i++)
                {
                    //log_info << "writing species position 2222  " <<i;
                    connectivityDataString += to_string(i);
                    connectivityDataString += "\n";
                }
                vertsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());
                //log_info << "writing connectivity 2222";
            }

            //offsets
            auto vertsNodeDataArrayOffsets = vertsNode.append_child("DataArray");
            vertsNodeDataArrayOffsets.append_attribute("type") = "Int32";
            vertsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
            vertsNodeDataArrayOffsets.append_attribute("format") = "ascii";

            if (config->pmpi->is_master())
            {
                string offsetsDataString = "\n";
                int iOffset = 0;
                for (int i = 0; i < particleNumberGlobal; i++)
                {
                    iOffset += 1;
                    offsetsDataString += to_string(iOffset);
                    offsetsDataString += "\n";
                }
                vertsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());
            }



            //Types
            auto vertsNodeDataArrayTypes = vertsNode.append_child("DataArray");
            vertsNodeDataArrayTypes.append_attribute("type") = "Int32";
            vertsNodeDataArrayTypes.append_attribute("Name") = "types";
            vertsNodeDataArrayTypes.append_attribute("format") = "ascii";

            if (config->pmpi->is_master())
            {
                string typesDataString = "\n";
                for (int i = 0; i < particleNumberGlobal; i++)
                {
                    typesDataString += "1";
                    typesDataString += "\n";
                }
                vertsNodeDataArrayTypes.text().set(typesDataString.c_str());
            }



            //get all velocity from all processes
            if (config->debugLevel > 1) log_info << "writing species vtk file: gather velocity";
            double* velocity = new double[particleNumberGlobal * 3];
            double* velocityTemp = new double[particleNumberLocal * 3];
            for (int i = 0; i < particleNumberLocal * 3; i++)
            {
                velocityTemp[i] = 0.0;
            }

            offset = new int[config->pmpi->get_size()];
            count = new int[config->pmpi->get_size()];
            offset[0] = 0;
            for (int iRank = 0; iRank < config->pmpi->get_size() - 1; iRank++)
            {
                offset[iRank + 1] = offset[iRank] + 3 * macroParticleNumber[i_species * config->pmpi->get_size() + iRank];
            }
            for (int iRank = 0; iRank < config->pmpi->get_size(); iRank++)
            {
                count[iRank] = 3 * macroParticleNumber[i_species * config->pmpi->get_size() + iRank];
            }


            for (int iRank = 0; iRank < config->pmpi->get_size(); iRank++)
            {
                int iParticle = 0;
                for (int i_volume = 0; i_volume < species->particles.size(); i_volume++)
                {
                    auto& cell_particles = species->particles[i_volume];
                    for (int i_particle = 0; i_particle < cell_particles.data.size(); i_particle++)
                    {
                        auto& particle = cell_particles.data[i_particle];
                        velocityTemp[iParticle * 3 + 0] = particle.velocity[0];
                        velocityTemp[iParticle * 3 + 1] = particle.velocity[1];
                        velocityTemp[iParticle * 3 + 2] = particle.velocity[2];

                        iParticle++;
                    }
                }
            }
            config->pmpi->gatherv_double(velocityTemp, 3*particleNumberLocal, velocity, count, offset, 0);

            //velocity
            auto cellDataNode = pieceNode.append_child("PointData");

            //log_info << "writing species velocity " << vtpFileName;
            if (config->debugLevel > 1) log_info << "writing species vtk file: write velocity";

            auto cellDataNodeNodeDataArrayVelocity = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayVelocity.append_attribute("Name") = "velocity";
            cellDataNodeNodeDataArrayVelocity.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayVelocity.append_attribute("NumberOfComponents") = "3";
            cellDataNodeNodeDataArrayVelocity.append_attribute("format") = "ascii";

            if (config->pmpi->is_master())
            {
                string velocityDataString = "\n";
                for (int iParticle = 0; iParticle < particleNumberGlobal; iParticle++)
                {
                    velocityDataString += to_string(velocity[iParticle * 3 + 0]);
                    velocityDataString += " ";
                    velocityDataString += to_string(velocity[iParticle * 3 + 1]);
                    velocityDataString += " ";
                    velocityDataString += to_string(velocity[iParticle * 3 + 2]);
                    velocityDataString += "\n";
                }
                cellDataNodeNodeDataArrayVelocity.text().set(velocityDataString.c_str());
            }
            delete[] velocityTemp;
            delete[] velocity;
            delete[] count;
            delete[] offset;


            if (config->pmpi->is_master())
            {
                if (config->debugLevel > 1) log_info << "writing species vtk file: save vtu file " << vtpFileName;

                doc.save_file(vtpFileName.c_str());
            }

            config->pmpi->barrier();

            log_info << "done writing species: " << species->species_param.species_name;

        }



    };

    void writeVtpBoundaryFieldGlobal(string timeString)
    {
        if (config->dimension == 3)
        {
            //writeVtpBoundaryFieldGlobal3d(timeString);
            writeVtuBoundaryFieldGlobal3d(timeString);
        }
        else if (config->dimension == 2)
        {
            //writeVtpBoundaryFieldGlobal2d(timeString);
            writeVtuBoundaryFieldGlobal2d(timeString);
        }
        else
        {
            log_error << "the dimension not support: " << config->dimension;
            exit(0);
        }
    }

    void writeVtpBoundaryFieldGlobal3d(string timeString)
    {
        EasyMesh::ParallelMesh* localMesh = static_cast<EasyMesh::ParallelMesh*>(mesh);
        EasyMesh::Mesh* globalMesh = localMesh->serialMesh;

        std::string vtpFileName = "result/vtk/boundaryField_" + timeString + ".vtp";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "PolyData";

        auto polyDataNode = vtkFileNode.append_child("PolyData");

        auto pieceNode = polyDataNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfVerts") = "0";
        pieceNode.append_attribute("NumberOfLines") = "0";
        pieceNode.append_attribute("NumberOfStrips") = "0";
        pieceNode.append_attribute("NumberOfPolys") = to_string(io->boundaryFields[0]->getFaceNumber()).c_str();

        //write points: 写网格点、面等的时候要用globalMesh
        if (config->debugLevel > 1) log_info << "writing point";
        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());

        //write face group polys(triangles): connectivity
        if (config->debugLevel > 1) log_info << "writing face";
        auto polysNode = pieceNode.append_child("Polys");

        auto polysNodeDataArrayConnectivity = polysNode.append_child("DataArray");
        polysNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        polysNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        polysNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;
            if (config->debugLevel > 1) log_info << "faceGroup name: " << faceGroupName;

            if (faceGroupName == "parallelFace")
            {
                continue;
            }

            for (int i = 0; i < faceGroup.faceList.size(); i++)
            {
                int iFace = faceGroup.faceList[i];
                EasyMesh::Face& face = globalMesh->faceList[iFace];
                connectivityDataString += to_string(face.pointList[0]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[1]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[2]);
                connectivityDataString += "\n";
            }
        }
        polysNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        //write face group polys(triangles): offset
        if (config->debugLevel > 1) log_info << "writing offset";
        auto polysNodeDataArrayOffsets = polysNode.append_child("DataArray");
        polysNodeDataArrayOffsets.append_attribute("type") = "Int32";
        polysNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        polysNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;

            if (faceGroupName == "parallelFace")
            {
                continue;
            }

            for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
            {
                iOffset += 3;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
        }
        polysNodeDataArrayOffsets.text().set(offsetsDataString.c_str());


        //PolyData存储的是面，以及面上的数据，因此这里的cell指的是Polys表示的面，值在面心
        auto cellDataNode = pieceNode.append_child("CellData");

        //将不同的BoundaryFIeldCollection写入不同的DataArray
        for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
        {
            if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
            auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayFlux.append_attribute("Name") = io->boundaryFields[iBoundaryField]->name.c_str();
            cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
            cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

            string fluxDataString = "\n";
            //write cell data
            if (config->debugLevel > 1) log_info << "globalMesh->faceGroupMap size: " << globalMesh->faceGroupMap.size();
            for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
            {
                auto faceGroupName = iter->first;

                if (faceGroupName == "parallelFace")
                {
                    continue;
                }

                if (faceGroupName == "")
                {
                    //临时跳过这个问题，后续要检查在哪里增加了空的faceGroup
                    log_error << "no faceGroupName";
                    continue;
                }
                if (config->debugLevel > 1) log_info << faceGroupName << " begin";

                auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(faceGroupName)->second;
                //boundaryField.set_value(0.0);
                for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                {
                    //Point& point = mesh->pointList[iPoint];
                    //potentialDataString += to_string(iPoint);
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << boundaryField.data_average_global[iFace];
                    ss >> str;

                    fluxDataString += str;
                    fluxDataString += "\n";
                }
                if (config->debugLevel > 1) log_info << faceGroupName << " end";

            }
            cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
        }
        if (config->debugLevel > 1) log_info << "end writing face data ";

        doc.save_file(vtpFileName.c_str());

    };

    void writeVtuBoundaryFieldGlobal3d(string timeString)
    {
        EasyMesh::ParallelMesh* localMesh = static_cast<EasyMesh::ParallelMesh*>(mesh);
        EasyMesh::Mesh* globalMesh = localMesh->serialMesh;

        std::string vtuFileName = "result/vtk/boundaryField_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto polyDataNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = polyDataNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(io->boundaryFields[0]->getFaceNumber()).c_str();

        //write points: 写网格点、面等的时候要用globalMesh
        if (config->debugLevel > 1) log_info << "writing point";
        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());

        //write face group polys(triangles): connectivity
        if (config->debugLevel > 1) log_info << "writing face";
        auto polysNode = pieceNode.append_child("Cells");

        auto polysNodeDataArrayConnectivity = polysNode.append_child("DataArray");
        polysNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        polysNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        polysNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;
            if (config->debugLevel > 1) log_info << "faceGroup name: " << faceGroupName;

            if (faceGroupName == "parallelFace")
            {
                continue;
            }

            for (int i = 0; i < faceGroup.faceList.size(); i++)
            {
                int iFace = faceGroup.faceList[i];
                EasyMesh::Face& face = globalMesh->faceList[iFace];
                connectivityDataString += to_string(face.pointList[0]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[1]);
                connectivityDataString += " ";
                connectivityDataString += to_string(face.pointList[2]);
                connectivityDataString += "\n";

                //log_info << "========================= iFace " << globalMesh->faceList.size() << " " << iFace;
                //log_info << "============= iPoint " << globalMesh->pointList.size() << " " << face.pointList[0];
                //log_info << "============= iPoint " << globalMesh->pointList.size() << " " << face.pointList[1];
                //log_info << "============= iPoint " << globalMesh->pointList.size() << " " << face.pointList[2];
                //log_info << "=============================================================================== ";
            }
        }
        polysNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        //write face group polys(triangles): offset
        if (config->debugLevel > 1) log_info << "writing offset";
        auto polysNodeDataArrayOffsets = polysNode.append_child("DataArray");
        polysNodeDataArrayOffsets.append_attribute("type") = "Int32";
        polysNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        polysNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;

            if (faceGroupName == "parallelFace")
            {
                continue;
            }

            for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
            {
                iOffset += 3;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
        }
        polysNodeDataArrayOffsets.text().set(offsetsDataString.c_str());

        //write face group polys(triangles): type
        if (config->debugLevel > 1) log_info << "writing type";
        auto polysNodeDataArrayTypes = polysNode.append_child("DataArray");
        polysNodeDataArrayTypes.append_attribute("type") = "Int32";
        polysNodeDataArrayTypes.append_attribute("Name") = "types";
        polysNodeDataArrayTypes.append_attribute("format") = "ascii";

        string typesDataString = "\n";
        for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
        {
            auto faceGroupName = iter->first;
            auto& faceGroup = iter->second;

            if (faceGroupName == "parallelFace")
            {
                continue;
            }

            for (int iFace = 0; iFace < faceGroup.faceList.size(); iFace++)
            {
                typesDataString += "5";
                typesDataString += "\n";
            }
        }
        polysNodeDataArrayTypes.text().set(typesDataString.c_str());


        //boundary field data is located at cell center
        auto cellDataNode = pieceNode.append_child("CellData");

        //write boundary field data: data_global
        if (config->debugLevel > 1)
        {
            for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
            {
                if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
                auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
                cellDataNodeNodeDataArrayFlux.append_attribute("Name") = (io->boundaryFields[iBoundaryField]->name + "_oneStep").c_str();
                cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
                cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
                cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

                string fluxDataString = "\n";
                //write cell data
                if (config->debugLevel > 1) log_info << "globalMesh->faceGroupMap size: " << globalMesh->faceGroupMap.size();
                for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
                {
                    auto faceGroupName = iter->first;

                    if (faceGroupName == "parallelFace")
                    {
                        continue;
                    }

                    if (faceGroupName == "")
                    {
                        //临时跳过这个问题，后续要检查在哪里增加了空的faceGroup
                        log_error << "no faceGroupName";
                        continue;
                    }
                    if (config->debugLevel > 1) log_info << faceGroupName << " begin";

                    auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(faceGroupName)->second;
                    //boundaryField.set_value(0.0);
                    for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                    {
                        //Point& point = mesh->pointList[iPoint];
                        //potentialDataString += to_string(iPoint);
                        std::stringstream ss;
                        std::string str;

                        ss.clear();
                        ss << std::setprecision(14) << boundaryField.data_global[iFace];
                        ss >> str;

                        fluxDataString += str;
                        fluxDataString += "\n";
                    }
                    if (config->debugLevel > 1) log_info << faceGroupName << " end";

                }
                cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
            }
        }
        

        //write boundary field data: data_average_global
        for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
        {
            if (io->boundaryFields[iBoundaryField]->name == "boundaryId")
            {
                continue;
            }

            if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
            auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayFlux.append_attribute("Name") = io->boundaryFields[iBoundaryField]->name.c_str();
            cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
            cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

            string fluxDataString = "\n";
            //write cell data
            if (config->debugLevel > 1) log_info << "globalMesh->faceGroupMap size: " << globalMesh->faceGroupMap.size();
            for (auto iter = globalMesh->faceGroupMap.begin(); iter != globalMesh->faceGroupMap.end(); iter++)
            {
                auto faceGroupName = iter->first;

                if (faceGroupName == "parallelFace")
                {
                    continue;
                }

                if (faceGroupName == "")
                {
                    //临时跳过这个问题，后续要检查在哪里增加了空的faceGroup
                    log_error << "no faceGroupName";
                    continue;
                }
                if (config->debugLevel > 1) log_info << faceGroupName << " begin";

                auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(faceGroupName)->second;
                //boundaryField.set_value(0.0);
                for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                {
                    //Point& point = mesh->pointList[iPoint];
                    //potentialDataString += to_string(iPoint);
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << boundaryField.data_average_global[iFace];
                    ss >> str;

                    fluxDataString += str;
                    fluxDataString += "\n";
                }
                if (config->debugLevel > 1) log_info << faceGroupName << " end";

            }
            cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
        }

        if (config->debugLevel > 1) log_info << "end writing face data ";

        doc.save_file(vtuFileName.c_str());

    };

    void writeVtpBoundaryFieldGlobal2d(string timeString)
    {
        EasyMesh::ParallelMesh* localMesh = static_cast<EasyMesh::ParallelMesh*>(mesh);
        EasyMesh::Mesh* globalMesh = localMesh->serialMesh;

        std::string vtpFileName = "result/vtk/boundaryField_" + timeString + ".vtp";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "PolyData";

        auto polyDataNode = vtkFileNode.append_child("PolyData");

        auto pieceNode = polyDataNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfVerts") = "0";
        pieceNode.append_attribute("NumberOfLines") = "0";
        pieceNode.append_attribute("NumberOfStrips") = "0";
        pieceNode.append_attribute("NumberOfPolys") = to_string(io->boundaryFields[0]->getFaceNumber()).c_str();

        //write points: 写网格点、面等的时候要用globalMesh
        if (config->debugLevel > 1) log_info << "writing point";
        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());

        //write face group polys(triangles): connectivity
        if (config->debugLevel > 1) log_info << "writing face";
        auto polysNode = pieceNode.append_child("Polys");

        auto polysNodeDataArrayConnectivity = polysNode.append_child("DataArray");
        polysNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        polysNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        polysNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
        {
            auto lineGroupName = iter->first;
            auto& lineGroup = iter->second;
            if (config->debugLevel > 1) log_info << "lineGroup name: " << lineGroupName;

            if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
            {
                continue;
            }

            for (int i = 0; i < lineGroup.lineList.size(); i++)
            {
                int iLine = lineGroup.lineList[i];
                EasyMesh::Line& line = globalMesh->lineList[iLine];
                connectivityDataString += to_string(line.pointList[0]);
                connectivityDataString += " ";
                connectivityDataString += to_string(line.pointList[1]);
                connectivityDataString += "\n";
            }
        }
        polysNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        //write line group polys(line): offset
        if (config->debugLevel > 1) log_info << "writing offset";
        auto polysNodeDataArrayOffsets = polysNode.append_child("DataArray");
        polysNodeDataArrayOffsets.append_attribute("type") = "Int32";
        polysNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        polysNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
        {
            auto lineGroupName = iter->first;
            auto& lineGroup = iter->second;

            if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
            {
                continue;
            }

            for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
            {
                iOffset += 2;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
        }
        polysNodeDataArrayOffsets.text().set(offsetsDataString.c_str());


        //PolyData存储的是面，以及面上的数据，因此这里的cell指的是Polys表示的面，值在面心
        auto cellDataNode = pieceNode.append_child("CellData");

        //将不同的BoundaryFIeldCollection写入不同的DataArray
        //只写出边界场的平均值
        for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
        {
            if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
            auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayFlux.append_attribute("Name") = io->boundaryFields[iBoundaryField]->name.c_str();
            cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
            cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

            string fluxDataString = "\n";
            if (config->debugLevel > 1) log_info << "globalMesh->lineGroupMap size: " << globalMesh->lineGroupMap.size();
            for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
            {
                auto lineGroupName = iter->first;

                if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
                {
                    continue;
                }

                if (lineGroupName == "")
                {
                    //临时跳过这个问题，后续要检查在哪里增加了空的lineGroup
                    log_error << "no lineGroupName: " << lineGroupName;
                    continue;
                }
                if (config->debugLevel > 1) log_info << lineGroupName << " begin";

                auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(lineGroupName)->second;
                //boundaryField.set_value(0.0);
                for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                {
                    //Point& point = mesh->pointList[iPoint];
                    //potentialDataString += to_string(iPoint);
                    fluxDataString += to_string(boundaryField.data_average_global[iFace]);
                    fluxDataString += "\n";
                }
                if (config->debugLevel > 1) log_info << lineGroupName << " end";

            }
            cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
        }
        if (config->debugLevel > 1) log_info << "end writing face data ";

        doc.save_file(vtpFileName.c_str());

    };

    void writeVtuBoundaryFieldGlobal2d(string timeString)
    {
        EasyMesh::ParallelMesh* localMesh = static_cast<EasyMesh::ParallelMesh*>(mesh);
        EasyMesh::Mesh* globalMesh = localMesh->serialMesh;

        std::string vtuFileName = "result/vtk/boundaryField_" + timeString + ".vtu";

        xml_document doc;
        auto vtkFileNode = doc.append_child("VTKFile");
        vtkFileNode.append_attribute("type") = "UnstructuredGrid";

        auto polyDataNode = vtkFileNode.append_child("UnstructuredGrid");

        auto pieceNode = polyDataNode.append_child("Piece");
        pieceNode.append_attribute("NumberOfPoints") = to_string(globalMesh->pointList.size()).c_str();
        pieceNode.append_attribute("NumberOfCells") = to_string(io->boundaryFields[0]->getFaceNumber()).c_str();

        //write points: 写网格点、面等的时候要用globalMesh
        if (config->debugLevel > 1) log_info << "writing point";
        auto pointsNode = pieceNode.append_child("Points");
        auto pointsNodeDataArray = pointsNode.append_child("DataArray");
        pointsNodeDataArray.append_attribute("type") = "Float64";
        pointsNodeDataArray.append_attribute("NumberOfComponents") = "3";
        pointsNodeDataArray.append_attribute("format") = "ascii";

        string pointDataString = "\n";
        for (int iPoint = 0; iPoint < globalMesh->pointList.size(); iPoint++)
        {
            EasyMesh::Point& point = globalMesh->pointList[iPoint];
            std::stringstream ss;
            std::string str;

            ss.clear();
            ss << std::setprecision(14) << point.position[0];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[1];
            ss >> str;
            pointDataString += str;
            pointDataString += " ";

            ss.clear();
            ss << std::setprecision(14) << point.position[2];
            ss >> str;
            pointDataString += str;
            pointDataString += "\n";
        }
        pointsNodeDataArray.text().set(pointDataString.c_str());

        //write cells: line
        if (config->debugLevel > 1) log_info << "writing face";
        auto cellsNode = pieceNode.append_child("Cells");

        auto cellsNodeDataArrayConnectivity = cellsNode.append_child("DataArray");
        cellsNodeDataArrayConnectivity.append_attribute("type") = "Int32";
        cellsNodeDataArrayConnectivity.append_attribute("Name") = "connectivity";
        cellsNodeDataArrayConnectivity.append_attribute("format") = "ascii";

        string connectivityDataString = "\n";
        for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
        {
            auto lineGroupName = iter->first;
            auto& lineGroup = iter->second;
            if (config->debugLevel > 1) log_info << "lineGroup name: " << lineGroupName;

            if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
            {
                continue;
            }

            for (int i = 0; i < lineGroup.lineList.size(); i++)
            {
                int iLine = lineGroup.lineList[i];
                EasyMesh::Line& line = globalMesh->lineList[iLine];
                connectivityDataString += to_string(line.pointList[0]);
                connectivityDataString += " ";
                connectivityDataString += to_string(line.pointList[1]);
                connectivityDataString += "\n";
            }
        }
        cellsNodeDataArrayConnectivity.text().set(connectivityDataString.c_str());

        //write line group polys(line): offset
        if (config->debugLevel > 1) log_info << "writing offset";
        auto cellsNodeDataArrayOffsets = cellsNode.append_child("DataArray");
        cellsNodeDataArrayOffsets.append_attribute("type") = "Int32";
        cellsNodeDataArrayOffsets.append_attribute("Name") = "offsets";
        cellsNodeDataArrayOffsets.append_attribute("format") = "ascii";

        string offsetsDataString = "\n";
        int iOffset = 0;
        for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
        {
            auto lineGroupName = iter->first;
            auto& lineGroup = iter->second;

            if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
            {
                continue;
            }

            for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
            {
                iOffset += 2;
                offsetsDataString += to_string(iOffset);
                offsetsDataString += "\n";
            }
        }
        cellsNodeDataArrayOffsets.text().set(offsetsDataString.c_str());


        //write : type
        if (config->debugLevel > 1) log_info << "writing type";
        auto cellsNodeDataArrayTypes = cellsNode.append_child("DataArray");
        cellsNodeDataArrayTypes.append_attribute("type") = "Int32";
        cellsNodeDataArrayTypes.append_attribute("Name") = "types";
        cellsNodeDataArrayTypes.append_attribute("format") = "ascii";

        string typesDataString = "\n";
        for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
        {
            auto lineGroupName = iter->first;
            auto& lineGroup = iter->second;

            if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
            {
                continue;
            }

            for (int iLine = 0; iLine < lineGroup.lineList.size(); iLine++)
            {
                typesDataString += "3";
                typesDataString += "\n";
            }
        }
        cellsNodeDataArrayTypes.text().set(typesDataString.c_str());



        //PolyData存储的是面，以及面上的数据，因此这里的cell指的是Polys表示的面，值在面心
        auto cellDataNode = pieceNode.append_child("CellData");

        //write boundary field data: data_global
        if (config->debugLevel > 1)
        {
            for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
            {
                if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
                auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
                cellDataNodeNodeDataArrayFlux.append_attribute("Name") = (io->boundaryFields[iBoundaryField]->name + "_oneStep").c_str();
                cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
                cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
                cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

                string fluxDataString = "\n";
                if (config->debugLevel > 1) log_info << "globalMesh->lineGroupMap size: " << globalMesh->lineGroupMap.size();
                for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
                {
                    auto lineGroupName = iter->first;

                    if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
                    {
                        continue;
                    }

                    if (lineGroupName == "")
                    {
                        //临时跳过这个问题，后续要检查在哪里增加了空的lineGroup
                        log_error << "no lineGroupName: " << lineGroupName;
                        continue;
                    }
                    if (config->debugLevel > 1) log_info << lineGroupName << " begin";

                    auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(lineGroupName)->second;
                    //boundaryField.set_value(0.0);
                    for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                    {
                        //Point& point = mesh->pointList[iPoint];
                        //potentialDataString += to_string(iPoint);
                        std::stringstream ss;
                        std::string str;

                        ss.clear();
                        ss << std::setprecision(14) << boundaryField.data_global[iFace];
                        ss >> str;

                        fluxDataString += str;
                        fluxDataString += "\n";
                    }
                    if (config->debugLevel > 1) log_info << lineGroupName << " end";

                }
                cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
            }
        }
        

        //write boundary field data: data_average_global
        for (int iBoundaryField = 0; iBoundaryField < io->boundaryFields.size(); iBoundaryField++)
        {
            if (io->boundaryFields[iBoundaryField]->name == "boundaryId")
            {
                continue;
            }

            if (config->debugLevel > 1) log_info << "writing face data: " << io->boundaryFields[iBoundaryField]->name;
            auto cellDataNodeNodeDataArrayFlux = cellDataNode.append_child("DataArray");
            cellDataNodeNodeDataArrayFlux.append_attribute("Name") = io->boundaryFields[iBoundaryField]->name.c_str();
            cellDataNodeNodeDataArrayFlux.append_attribute("type") = "Float64";
            cellDataNodeNodeDataArrayFlux.append_attribute("NumberOfComponents") = "1";
            cellDataNodeNodeDataArrayFlux.append_attribute("format") = "ascii";

            string fluxDataString = "\n";
            if (config->debugLevel > 1) log_info << "globalMesh->lineGroupMap size: " << globalMesh->lineGroupMap.size();
            for (auto iter = globalMesh->lineGroupMap.begin(); iter != globalMesh->lineGroupMap.end(); iter++)
            {
                auto lineGroupName = iter->first;

                if (lineGroupName == "parallelFace" || lineGroupName == "parallelLine")
                {
                    continue;
                }

                if (lineGroupName == "")
                {
                    //临时跳过这个问题，后续要检查在哪里增加了空的lineGroup
                    log_error << "no lineGroupName: " << lineGroupName;
                    continue;
                }
                if (config->debugLevel > 1) log_info << lineGroupName << " begin";

                auto& boundaryField = io->boundaryFields[iBoundaryField]->boundaryFieldMap.find(lineGroupName)->second;
                //boundaryField.set_value(0.0);
                for (int iFace = 0; iFace < boundaryField.globalSize(); iFace++)
                {
                    //Point& point = mesh->pointList[iPoint];
                    //potentialDataString += to_string(iPoint);
                    std::stringstream ss;
                    std::string str;

                    ss.clear();
                    ss << std::setprecision(14) << boundaryField.data_average_global[iFace];
                    ss >> str;

                    fluxDataString += str;
                    fluxDataString += "\n";
                }
                if (config->debugLevel > 1) log_info << lineGroupName << " end";

            }
            cellDataNodeNodeDataArrayFlux.text().set(fluxDataString.c_str());
        }


        if (config->debugLevel > 1) log_info << "end writing face data ";

        doc.save_file(vtuFileName.c_str());

    };

private:
    EasyMesh::ParallelMesh* mesh;

    EasyMesh::Mesh* globalMesh;

    Config* config;
    IO* io;
};


#endif
